Phase diagram of a two-dimensional system with anomalous liquid properties 
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Using Monte Carlo simulation techniques, we calculate the phase diagram for a square shoulder- 
square well potential in two dimensions that has been previously shown to exhibit liquid anomalies 
consistent with a metastable liquid-liquid critical point. We consider the liquid, gas and five crystal 
phases, and find that all the melting lines are first order, despite a small range of metastability. One 
melting line exhibits a temperature maximum, as well as a pressure maximum that implies inverse 
melting over a small range in pressure. 



I. INTRODUCTION 



Core-softened potentials were first used by Stell, Hem- 
mer and coworkers in a lattice gas system to discuss the 
isostructural solid-solidphase transition that ends in a 
second critical point [l|-y| • Core-softened potentials were 
also used to study single-component systems in a liq- 
uid state, such as liquid metals |4J— fill] . They have been 
also used to study liquid anomalies in ID |12h14J and 
2D [lMl|. Calculations in Ref. 0, HJ show that a 
core-softened potential can be considered as a realistic 
first-order approximation for the real interaction between 
water molecules resulting from averaging over the angu- 
lar part. 

Interest in the study of liquid-liquid (L-L) phase tran- 
sitions in single component systems grew dramatically 
after such a transition and accompanying critical point 
were proposed for water as an explanation for its anoma- 
lous properties [21| . Various studies have been done to 



understand the L-L phase transition and associated phe- 
nomena. Some of these studies focus on the "two-liquid" 
model to explain liquid properties |22h24| . Other stud- 
ies were based on using anisotropic potentials [25|, |26| . 
Franzese et al. showed that the liquid-liquid phase tran- 
sition and accompanying critical point can also arise from 
an isotropic interaction potential with two characteristic 
distances (hard-core and soft-core) [27j |. In this work, 
the authors reported in 3D molecular dynamics (MD) 
simulations the existence of two liquid phases, the low- 
density liquid phase and the high-density liquid phase, 
and showed that these two phases can occur in the sys- 
tem with no density anomaly. On the other hand, 2D 
simulation studies re pro duce the density anomaly but no 
second critical point [13J . For a review of unusual behav- 
ior of isotropic potentials with two energy scales in 2D, 
see Ref. [H|. 

Scala and coworkers [18[ carried out MD simulations 
in 2D of the square-shoulder square- well (SSSW) poten- 
tial shown in Fig. Q] to study liquid anomalies. Buldyrev 
et al [H continued with the SSSW model in 2D and 3D 
in order to study liquid-liquid phase transitions. For the 
2D system, they produced a phase diagram showing liq- 
uid anomalies in relation to approximate crystallization 
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FIG. 1: The pair potential used in this study is an isotropic 
step potential with hard-core diameter a. b — y/~2a is the soft- 
core distance, and c = \/Ha is the attractive distance limit, r 
is the distance between two particles and e is the bond energy. 



lines for a range in pressure P and temperature T near a 
potential L-L critical point. Their phase diagram shows 
the gas-liquid coexistence curve and crystallization lines 
for a low density triangular and higher density square 
crystal. It also shows the first critical point and the hy- 
pothetical position of the second critical point, which co- 
incides with the crossing of the two crystallization lines. 
Thus, unavoidable crystallization renders the L-L crit- 
ical point not directly observable, or obscured. Their 
crystallization lines were determined from examining the 
behavior of the pressure, structure and dynamics along 
isochores. They are estimates of the limit of liquid sta- 
bility, or rather the limit of metastability, with respect to 
the crystal, rather than thermodynamically determined 
coexistence lines. As the system is two-dimensional, the 
nature of the crystallization transition is also under ques- 
tion, in so far that in two dimensions, crystallization can 
proceed in a continuous way via a hexatic phase rather 
than through a first-order phase transition. 

In the present work, we carry out free energy calcula- 
tions, based primarily on Monte Carlo (MC) simulation, 
to determine the coexistence conditions between the liq- 



uid and crystal phases for a wide range of P and T, in- 
cluding the smaller range presented in Ref . [29| . In doing 
so, we find two low density crystal phases not previously 
reported for the model. We find that all the transitions 
are at least weakly first-order. The crystallization lines 
reported in Ref. [29| are below our calculated melting 
lines. Additionally, the square crystal shows a maximum 
temperature in its melting curve, as well as a maximum 
in pressure. Thus, the present model is a useful one 
for studying the rare phenomenon of inverse melting, in 
which the liquid may freeze to the crystal upon heating. 
This paper is organized as follows. In Section II, we 
discuss all the free energy and computer simulation tech- 
niques used in carrying out this work. In section III, we 
show our results. In Section IV we present a discussion 
and we give our conclusions in Section V. 



II. METHODS 

A. Model and simulations 

The model we study is the step pair potential shown 
in Fig. [T] As we are carrying out our studies in two 
dimensions, the model describes disks with a hard-core 
diameter a and an attractive well extending out to a 
radial distance of c = ^/3a. The attractive well itself 
contains a shoulder, with a pair interaction energy of 
— e/2 for a < r < b and energy of — e for b < r < c. 
The parameter b was originally chosen to be v2er so that 
there would exist a low density triangular (LDT) phase 
and a higher density square (S) phase with the same po- 
tential energy per particle of — 3e, i.e. two energetically 
degenerate phases of well separated densities [29j ■ The 
idea behind this was to allow for distinct liquid states, 
one based on square packing and the other on the more 
open triangular packing, in analogy to what is thought 
to be the case for water. At high pressure the system 
ultimately must form the close-packed triangular phase 
(HDT), with potential energy per particle of — 1.5e. We 
find two additional crystals, phases A and Z, with per 
particle energies — 3.25e and — 3.5e, respectively. The 
various crystal phases are depicted in Fig. [2j Our goal 
is to calculate coexistence lines between the five crystal 
phases, the liquid (L) and the gas (G). 

The liquid-state properties of the model were exten- 
sively studied in Ref. [29| using discrete MD simula- 
tion. The S and LDT crystallization lines were deter- 
mined in that work from pressure isochores and from 
direct observation of crystal-like structural and dynami- 
cal behavior. Here, we calculate the crystal coexistence 
lines using free energy techniques that employ for the 
most part MC simulations performed at constant par- 
ticle number N, P and T, i.e., in the NPT ensem- 
ble [20]. Depending on the phase, the pressure is kept 
constant by changing the volume isotropically (for L, S, 
HDT and LDT), by allowing rectangular dimensions of 
the simulation cell to change length independently while 



maintaining a right angle (for A and Z), or by allowing 
the angle to change as well (as a check for all phases). 
The system sizes and box shapes are as follows: (L) 
N = 1020 and 986, square box; (S) N = 1024 square 
box, and N = 992 with rectangular box L y = 32L X /31; 
(HDT and LDT) N = 986, L y = 17y/3L x /29; (A) 
N = 952, L y = 28(sinl2° + l)Z x /(34 cos 12°) initially; 
(Z) N — 968, square box initially. The different box 
shapes (and hence number of particles) are used as con- 
sistency checks, and indeed we do not detect any differ- 
ence in the results based on the particular choice used. 
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FIG. 2: Illustration of the phases modelled: (a) the liquid 



(L), here shown as a small portion of a simulation in which 



distinct local packing environments are visible, (b) the square 
crystal (S), |(c)| the low-density triangular crystal (LDT), [(d)] 
the high-density triangular crystal (HDT), (e) the A crystal 
and (f) the Z crystal. Line segments for the crystal phases 
indicate a bond with energy —e and a dashed line segment 
one with energy —e/2. 



B. Solid-liquid and solid-solid coexistence 

First-order transition lines can be determined us- 
ing a method developed by Kofke to trace coexistence 
curves [3l|, [32[. Kofke refers to his method as Gibbs- 
Duhem integration, and it is based on the Clapeyron 
equation which describes the temperature dependence of 



the pressure at which two phases coexist, 



dP_ 
~dT 



As 
Av 



Ah 
TAJ? 



(1) 



where As is the molar entropy difference, Ah is the mo- 
lar enthalpy difference and Aw is the molar volume dif- 
ference between the two coexisting phases. Tracing the 
coexistence curve requires that one point on the coexis- 
tence curve be known and then the rest of the curve can 
be found by integration of Eq. [TJ in particular using the 
enthalpy since it is much easier to calculate than the en- 
tropy. We carry out the integration using a second-order 
predictor-corrector method [33|, |34| . 

To obtain the first coexistence point between the liq- 
uid and the S crystal, we first determine the respective 
equations of state along an isotherm by carrying out sev- 
eral NPT simulations. We choose fcsT/e = 0.55 so that 
we are above the L-G critical temperature, where fcs is 
the Boltzmann constant. Once the equations of state are 
known, we calculate the chemical potential \x for each 
phase as a function of number density p by integrating 
the pressure via [30, |35[ , 



0l*(p)=Pf(p*)+P 



~ d P 



P * P' 



£, (2) 

P 



where (3 = (ksT)' 1 and / is the Helmholtz free energy 
per particle calculated at a reference number density p* . 
To carry out the integration, we fit the liquid isotherm 
to Eq. [3] and the solid isotherm to Eg . HI [36l [37lj . 
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1 - a L p 



,(3) 
(4) 



where a;, s , 6/, s , and Q )S are the fit parameters. Integra- 
tion of Eq. [3J from zero to a density of interest yields the 
chemical potential of liquid, as given in Eq. [5] Similarly, 
integration of Eq. @] from a reference density to the den- 
sity of interest yields the chemical potential of solid, as 
given in Eq. © J2l III , 



PPi(p) = In 



pk 2 



bi/ai - ci/af + 1 



cip 



1 -cup 

ci/2a 2 + bip 

(l-a lP f ^ {1-aipf 

- (6,/oj - ci/2o? + 1), 

0li„(j>) = 2a s p + b s [\n(p) + 1} 

- [a s p*+b s ln{p*)-c s /p*] 
+ (3r(p*) + \n(A 2 p*)-l, 



dip 



(5) 



(6) 



an isotherm). / ox (p*) is the excess Helmholtz free energy 
per particle calculated at p* . 

For the liquid, Eq. [3J provides a good fit only up to 
p f=a 0.1, and so from p = to p\ = 0.09418 we use 
Eq. [5j and then integrate Eq. [2] numerically, using dif- 
ferent interpolation orders to estimate uncertainty. The 
equations of state for the liquid and the S crystal are 
shown in Fig. [3J 
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FIG. 3: Equations of state of the liquid (circles) and S crystal 
(diamonds) at feT/e = 0.55. The curves show fits according 
to Eqs. [3J and [1 with a t = -0.457818, b t = -9.18372, and 
ci = 33.4007 for the liquid (inset) and a s = 479.035, b s = 
-686.583, and c s = 246.067 for the crystal. 

We calculate the crystal reference Helmholtz free en- 
ergy using the Frenkel-Ladd method Q. In this method, 
a harmonic potential is added to the original system to 
define a new system potential energy, 



U x 



N 

U(r N ) + \Y,(n 

i=l 



ro,if 



(7) 



where fj is the position of particle i and r*o,i is its ideal 
lattice position, and U(r N ) is the unaltered system po- 
tential energy. U\ is such that at coupling parameter 
A = the original model is recovered and for sufficiently 
large A, the system behaves as an ideal Einstein crystal. 
A thermodynamic integration at a particular T and p is 
carried along A to determine the Helmholtz free energy 
difference between the Einstein crystal and the original 
model. The excess free energy per particle for the model 
is then expressed as |30 |. 
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Phu 
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CM 
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Hp*) 
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2N 



ln(AT) 
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where A = hj ' yj (2-KmkBT) is the de Broglie thermal 
wavelength, where it is assumed to equal unity since it 
plays no rule in locating the coexistence pressure (along 



where d — 2 is the dimensionality of the system, m = 
1 is the mass of the particle. The first term in Eq. [5] 
represents the free energy of the ideal (non-interacting) 



Einstein crystal, which is equal to, 



N 



In 



2/3 V/3A 



(9) 



where U(fQ ) is the potential energy of the crystal when 
all the atoms are at their ideal lattice positions, and A max 
is chosen such that, for A larger than A max , the mean- 
squared displacement (5r 2 )\ = ((rj — r 0yi ) 2 )\, where 
(. . . )\ indicates an ensemble average, for a system with 
fixed center of mass follows the following analytical ex- 
pression, 



(Si 



/Eins,A — 



N- 1 1 

N J\' 



(10) 



The second term in Eq. [5] represents the free energy dif- 
ference between the solid and the Einstein crystal, and 
can be calculated by integrating the mean-squared dis- 
placement obtained from simulations carried out with a 
fixed center of mass as follows [30, 39], 



AF 



CM 



N 



(6r 2 ) x dX. 



(11) 



This integration can be understood as gradually switch- 
ing on the coupling parameter to transform the solid into 
an Einstein crystal. For better accuracy, this integral can 
be transformed to 1301 , 



AF 



CM 



N 



ln(A„ 
ln(c) 



«+c 



#n(A + c)](A + c)(r 2 ) A , 



(121 



where c is a constant chosen to be 1 in this work. The 
integrand is shown in Fig. [4[ along with the curve for 
the ideal solid. We choose In (A max + 1) = 6.909, check- 
ing that using higher values yields no appreciable change 
in the final result. The integration is carried out using 
interpolations of different order in order to estimate un- 
certainty. 

The third, fourth and fifth terms in Eq. [5] correspond 
to the difference between the constrained (fixed center 
of mass) and unconstrained (non-fixed center of mass) 
solids. The last term in Eq. [5] is the free energy of the 
ideal gas per particle, which is given by, 



/3/ id = ln(p) - 1 + 



ln(27r7V) 

2N 



(13) 



Once the chemical potentials of the two phases are 
known, the coexistence point can be obtained from the 
intersection of the two chemical potential curves [3a, [37| . 

fj-i(p) and Hs(p) are used together with the equations 
of state to plot the chemical potentials of the two phases 
as functions of pressure, as we do in Fig. [SJ It is im- 
mediately apparent that fi(P) has nearly the same slope 
for both phases, and hence the location of the crossing 
is sensitive to errors in the various calculated quantities 
used to determine the curves. We note that the equa- 
tions of state are determined only to the point where the 
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FIG. 4: The mean-squared displacement transformed by 
Eq. [12] as a function of coupling parameter A calculated 
by computer simulation (solid curve is a guide to the eye). 
Dashed line is the theoretical value given by Eq. 1101 



metastable phase does not easily transform to the other 
phase. It is somewhat surprising that at the P for which 
either phase becomes unstable, Pa 2 /e ~ 3.49 for S and 
Pa 2 /e ~ 4.09 for L, the difference in chemical potential 
is very small, on the order of |/3A/z| ~ 0.01. 

As a check on the L-S coexistence conditions at 
ksT /e = 0.55, we perform an NVT (canonical ensem- 
ble) simulation with 10,000 particles initially placed on 
a square lattice with pa 2 = 0.786567, the p at which 
the system is expected to phase separate into L and S 
with equal numbers of particles in each phase, based on 
liquid and S coexistence densities of pi = 0.7677 and 
p x = 0.8064, respectively. Fig. [?5] shows a snapshot after 
running for 2 x 10 7 MC steps per particle, with dark sym- 
bols identifying particles belonging to the S phase [401 - 
|43| ] . Averaging over the last 5 x 10 6 MC steps per particle, 
the fraction of particles belonging to the S phase is 0.51. 

The above procedure is repeated (at lower T) for the 
other crystal phases to determine crystal-crystal coexis- 
tence lines. For two crystals, the slopes of n(P) are gen- 
erally quite different, which makes it easier to pinpoint 
the coexistence P. Similarly, at T less than the L-G criti- 
cal temperature, the procedure is repeated to find crystal 
sublimation lines after determining the equation of state 
for the gas. 

For the L-LDT melting line, we must additionally per- 
form an integration of the enthalpy H to lower T at a 
P above the critical pressure in order to avoid the L- 
G critical point. Specifically, we first integrate the liq- 
uid equation of state at ksT/e = 0.70 using Eq. [2] to 
Pa 2 /e = 0.05, and then calculate p(T) via [39|, 



p(T 2 ,P) _ fi(T 1 ,P) 
ksT2 ksTi 



'A 



H(T) 
Nk B T 2 



dT 7 



(14) 



noting that here, the T dependence of A must be taken 
into account. Equivalently, this amounts to using the po- 
tential energy instead of the thermal energy in calculating 



8 


1 i ' i 

-(a) 


1 i ' i ', 


6 




— Liquid phase 

- - Square crystal 


y'' 


4 






,- 










s 


1 i ' ' ' \ > 




CO. 2 




>'' 2.4 

=1 


- / \ 


- 






CO. 


s 









2 


/. 1 , . , 1 . 


■ 


-2 




3.6 , 4 


- 


.1 


'.i,i 


Po7e 

i.i. 




FIG. 5: Determination of a coexistence P between the L 
and S phases at ksT/e = 0.55. Panel (a) shows the chemical 
potential isotherms for the liquid (solid curve) and the square 
crystal (dashed curve). Inset shows a close-up of the crossing. 
In panel (b) we show the difference in chemical potential Ap 
between the two phases over the entire range of P for which 
the equations of state overlap, with dashed lines indicating 
upper and lower uncertainty estimates. 



H. For the LDT crystal, the reference free energy is cal- 
culated at Pa 2 /e — 0.05 after determining the density at 
that pressure to be pa 2 = 0.4780 ± 0.0015. In Fig. [7J we 
show H(T) for L and S as well as the resulting difference 
in /i between the phases. We repeat the calculation us- 
ing the the liquid equation of state at ksT/e = 0.55 as 
a check. Using the same procedure at ksT/e = 0.70, we 
carry out an evaluation of the melting temperature of the 
S phase at Pcr 2 /e = 0.15 [Fig.GJc)] and Pa 2 /e = 7.00 as 
a check on the accuracy of the coexistence line. 



C. L-HDT coexistence 

Using the Gibbs-Duhem integration method is not nec- 
essary (or possible) for tracing out the L-HDT melting 
line at high T, as over a certain range in P the system 




FIG. 6: Snapshot configuration obtained from an NVT simu- 
lation for 10,000 particles at k B T/e = 0.55 and p = 0.786567. 
Black symbols represent particles belonging to the S phase, 
while grey symbols represent the L phase. 



can fairly easily sample both states. Thus, to determine 
the coexistence P along an isotherm, we first locate a 
pressure Pq for which we can sample both states with 
reasonable statistics, as shown in Fig. |8l and determine 
the conditional Gibbs free energy from a histogram of the 
densities sampled during an NPT simulation, 



pAG(T ) P ;p) = -ln[P r (p)} ) 



(15) 



where P r (p) is the probability density of observing the 
system at a particular p. Here, we do not normalize our 
histograms as the normalization merely adds an incon- 
sequential shift. Pq already provides an estimate of the 
location of the coexistence pressure. The conditional free 
energy shown in Fig. [5] (black curve) exhibits a global 
minimum at high density (HDT) and a metastable one 
at low density (liquid). The free energy barrier between 
the two states is characteristic of a first order transi- 
tion. To more precisely locate the coexistence pressure, 
we reweight the histogram by applying a pressure shift, 



pAG(T,P';p) = f3AG(T,P ;p) 



N/3AP 
P 



(16) 



where c is a constant related to normalization and AP is 
the pressure shift that brings the two minima to the same 
level, as in Fig. [9] (red curve). The coexistence pressure 
is then equal to P' — Pq + AP. In practice, the shift 
we obtain is hardly perceptible on the scale of our plots, 
e.g., for the k B T/e = 5.0 case in Fig. 1 P a 2 /e = 50.0 
and AP<r 2 /e = -0.135, and for k B T/e = 1.0, P a a 2 /e = 
14.350 and APa 2 /e = -0.004. We note that the barrier 
does grow with decreasing T, and below fc^T/e w 0.5, 
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FIG. 7: (a) Enthalpy per particle for the liquid (circles) and 
LDT crystal (diamonds) along Pa 2 /e — 0.05. Here we have 
subtracted the ideal gas contribution to the energy, (b) The 
corresponding chemical potential difference between the L and 
LDT phases for the entire range in T of metastability. (c) The 
chemical potential difference between the L and S phases at 
Pa 2 /e = 0.15. 



both phases can stably exist for sufficiently long times 
in order to perform Gibbs-Duhem integration. Indeed 
below this T, it is not feasible to continue with histogram 
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FIG. 8: Sample time series of the number density near the 
L-HDT coexistence curve, with N — 986, ksT/e = 5.0 and 
Po = 50.0e/<7 2 . The system samples both the (lower density) 
liquid and the HDT crystal. 
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FIG. 9: Conditional Gibbs free energy as a function of p. 
At ksT/e = 5.0 and a pressure Po = 50.0e/<r 2 slightly above 
coexistence (solid curve, circles), the high density basin (HDT 
crystal) has a lower free energy than the low density (liquid) 
basin. Through Eq. 1161 an appropriate shift in the pressure 
locates the coexistence pressure, i.e., transforms the Po curve 
so that the liquid and HDT minima are at the same level to 
within precision of the data (dashed line, squares). 



reweighting without using some biasing potential within 
the MC simulations. 



D. G-L coexistence 

The G-L coexistence line can be determined by us- 
ing the Gibbs ensemble MC method developed by Pana- 
giotopoulos J44J. The Gibbs ensemble employs two sepa- 
rated subsystems (without the presence of an interface), 
where the total number of particles is fixed and the to- 
tal volume (in this case, area) of the two subsystems is 



also fixed; the total system as a whole evolves accord- 
ing to the canonical ensemble. The thermodynamic re- 
quirements for phase coexistence are that the tempera- 
ture, pressure, and chemical potential of the two coexist- 
ing phases must be equal and these requirements can be 
achieved by performing three different kinds of trial MC 
moves. First, particle displacement within each subsys- 
tem, second, volume fluctuations of the two subsystems, 
and third, transferring particles between the two subsys- 
tems. The advantage of using the Gibbs ensemble is that 
the system finds the densities of the coexistence phases 
without computing either the pressure or the chemical 
potential. 

Having obtained the coexistence densities at a series of 
T, the corresponding coexistence pressures can be esti- 
mated by applyin g th e virtual volume change method of 
Haresmiadis et al [451 ] . In this method, we perform sepa- 
rate NVT MC simulations of both the liquid and the gas 
at their respective coexistence densities (at a given T), 
and obtain the pressure via, 



k B T 

AV 



In 



y) cxp(-/3AC/)\ 



(17) 



where AU is the potential energy difference between 
a configuration with particle coordinates isotropically 
rescaled to accommodate a smaller virtual area V and 
the unaltered configuration with original area V, where 
V = V - AV and AV = O.ler 2 . Both phases give the 
same pressure to within error. 

However, as the temperature approaches the critical 
temperature Tc, G-L coexistence can no longer be dis- 
cerned in the Gibbs ensemble simulation. Our data for 
the G-L coexistence curve from the Gibbs ensemble ex- 
tend only to k B T/e = 0.50. Beyond this T, we extrap- 
olate according to the following procedure. We estimate 
Tc by fitting the density difference of the two coexisting 
phases to a scaling law [30, [36|, |46| , 



Pl -p g =A\T-T c f 



(18) 



where /? c is the critical exponent, which is equal to 0.125 
for a two-dimensional system, and A is a constant de- 
termined from the fit. To estimate the critical density 
PC-, we fit our results to the law of rectilinear diame- 
ters in mEi, 



Pl +Pg 



p c + B\T-T c \, 



(19) 



where B is a constant determined in the fit, and Tc is 
used from the fit in Eq. [TSJ The critical pressure Pc 
is estimated by fitting the vapor pressure curve to the 
Clausius-Clapeyron equation [46j . 
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c + f, 



(20) 



where C and D are constants determined in the fit. Pc is 
then calculated by substituting Tc obtained from Eq. [TH] 



in Eq.l20l From the fits, we obtain p c a 2 = 0.263 ±0.002, 
k B T c /e = 0.533 ± 0.002 and P C cr 2 /e = 0.019 ± 0.001. 
The uncertainties quoted here are based on uncertainties 
in the fit parameters and do not reflect any systematic 
error associated with the fact that we are extrapolating 
above k B T/e = 0.50, the highest T at which we have 
reliable Gibbs ensemble data. 



III. RESULTS 

Having assembled all of the individual coexistence 
curves, we present the phase diagram in the P-T plane in 
Fig. [10] and in the p-T plane in Fig. [TTJ The three pan- 
els of Fig. [TU]show progressively smaller ranges of P. In 
Fig. |10(b)[ dashed lines indicate metastable extensions of 
coexistence lines into the gas stability field (i.e., showing 
the phase diagram in the absence of the gas). As an aid 
to interpreting Fig. 111! we recall that under conditions 
of constant volume, the thermodynamic ground state is 
not necessarily a single phase, but is generally composed 
of two coexisting phases. The bottom panel of Fig. [TTJ 
shows the phase diagram in the absence of the gas phase. 



Fig. 10(a) shows a prominent S-L melting line temper- 
ature maximum at Per 2 /e = 5.24 ± 0.05 and k B T mllx /e = 
0.655 ± 0.005. At this point, according to Eq. [T] the 
molar volumes of the S crystal and liquid are equal. At 
higher P, the melt is more dense than the crystal, as in 
the familiar case of water and hexagonal ice. 

An even more exotic feature of the S-L melting line 
is the pressure maximum occurring near the HDT-S-L 
triple point at P max cr 2 / e = 7.98 ± 0.08 and k B T/e = 
0.450 ± 0.003. A close-up of this feature is shown in 
Fig. [TSJ At this point, according to Eq. [TJ the entropy 
of the S crystal and the liquid are equal, and for lower 
T along the curve, the melt has a lower entropy than 
the crystal. The presence of the pressure maximum in 
the melting curve allows for "inverse melting" [47| in the 
narrow range of P between the triple point and the max- 
imum, i.e. isobaric heating of the liquid results in crys- 
tallization. 

Given the numerical uncertainties in determining co- 
existence conditions and tracing out coexistence lines, we 
carry out a rough check by preforming three sets of simu- 
lations in the vicinity of the HDT-S-L triple point. Each 
set is a grid of 121 simulations for state points marked in 
Fig.[l2j For one set, the particles are initially positioned 
on the S lattice; for the second set, points are initially on 
the HDT lattice; high T liquid state configurations seed 
the third set of simulations. We run each simulation for 
5 x 10 7 MC steps per particle, and then indicate with the 
appropriate symbol in Fig.[T2]thc phase which the system 
spontaneously adopts. Potentially, since there are three 
simulations per state point, three symbols may appear, 
indicating stability or metastability of all three phases. 
Near the triple point, the simulations retain the starting 
phase, as expected, while deep within a stability field, all 
sets transform to the same phase. In this way, we crudely 
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FIG. 10: Phase diagram of the 2D model in the P-T plane, 
showing the liquid (L), gas (G) and crystal phases HDT, S, 
LDT, A and Z (see Fig. [2}. The panels show portions of the 
phase diagram at |(a)| high, |(b)| medium and |(c)| low P. The 
liquid-gas coexistence line terminates at a critical point at 
k B T c /e = 0.533 and P c a 2 /e = 0.0185 (filled circle). Dashed 
lines in (b) are metastable coexistence lines assuming the ab- 
sence of the gas phase. Initial coexistence points, i.e., start- 
ing points for Gibbs-Duhem integration, are indicated by cir- 
cles, while x 's show repeated coexistence calculations done as 
checks on the Gibbs-Duhem integration. 




FIG. 11: Phase diagram in the p-T plane, (a) The points 
along the G-L coexistence lines indicate results from Gibbs 
Ensemble simulations and the large filled orange circle shows 
our estimate of the G-L critical point based on an extrapola- 
tion described in the text. Panel (b) shows the phase diagram 
in the absence of the gas phase. The filled black circle shows 
the location of the obscured L-L critical point discussed in 
Ref. [H. 



map out the extent of metastability. 

It is difficult to directly confirm inverse melting on typ- 
ical simulation time scales, as the metastable phase is 
never far from the coexistence line. We aim to address 
this in future work. However, and while this is not a 
definitive check on the existence of inverse melting, the 
tendency for points exhibiting liquid metastability within 
the S stability field to track the curvature of the S-L melt- 
ing line is supportive of the existence of this phenomenon 
in the system, i.e., the lowest P point for each T for which 
the x and o simultaneously occur roughly form a curve 
with a maximum in P that tracks the shape of the S-L 
coexistence line. 

At lower P, we confirm the negative slope of the LDT- 
L melting line as reported already in Ref. 29] . Below the 
LDT-S-L triple point, we find that the new crystal phases 
A and Z both have reasonably large stability fields, as 
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FIG. 12: P-T phase diagram obtained at high pressure, near 
the L-S-HDT triple point. The grid ol points is obtained 
from three sets of simulations. Each of L, S and HDT is 
used to initialize a simulation set with N = 986, 992 and 
986, respectively. The final phase adopted from each set at 
each state point is indicated by a symbol: S, square; HDT, 
triangle; L, x. E.g., at low P and high T, both L and HDT 
transform to S, while near the triple point, each phase retains 
metastability. 



shown in Fig. |10(b)| and that the LDT crystal, having 
the lowest density of the crystal phases studied, occupies 
a rather small portion of the phase diagram. The A-S 
transition line is also negatively sloped, which together 
with the fact that the A phase has a lower density than 
S (see Fig.!!!]), implies through Eq. [T]that the entropy of 
S is larger than that of A. Indeed, the bonding distances 
required to form A are rather restrictive compared to the 
geometry of S, and this is reflected in the smaller range in 
p for which A is the single stable phase (again, compared 
to S). A similar argument holds when comparing Z to A. 

In Ref. [29(, the authors locate lines in the P-T plane 
that demarcate a limit to observing the liquid, i.e., where 
crystallization is practically unavoidable. Although their 
investigation into this aspect of the model was not ex- 
haustive, the character of crystallization was possibly 
suggestive of continuous crystallization seen in other two 
dimensional systems. We plot these lines within the 
appropriate portion of our calculated phase diagram in 
Fig. [T3l We see that the crystallization lines occur well 
below our calculated first-order melting lines, and there- 
fore occur at conditions for which there is a gap in crystal 
and liquid chemical potential. However, this does merit a 
closer look at the crystallization process, especially near 
the apparent limit of liquid metastability. Also in Fig.lT51 
we plot the location of what might be termed the ob- 
scured L-L critical point at low T that appears to be 
responsible for the liquid anomalies reported in Ref. |29j , 
but which is unobservable owing to unavoidable nucle- 
ation. Within uncertainty, this obscured critical point 
falls on the S-LDT coexistence line. 

We estimate the density of the obscured critical point 
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FIG. 13: Comparison of our phase diagram with previ- 
ously reported system properties. Red curves are taken from 
Ref. [2y] and represent crystallization lines (+), locus of tem- 
peratures of maximum density along isobars (TMD, circles), 
pressures of maximum diffusivity along isotherms (D max , 
diamonds), maxima of isothermal compressibility (^Tmax, 
squares) and G-L coexistence. Also shown are the G-L criti- 
cal point (filled red circle) and the obscured L-L critical point 
(large hashed circle). Dotted lines show pressure along iso- 
chores. All other symbols as in Fig. 1101 We note that we 
determine the location of the G-L critical point from an ex- 
trapolation of data above fc_eT/e = 0.50, while the one re- 
ported in Ref. [20] is based on inflection points of pressure 
isotherms. The previously reported crystallization lines fall 
within the presently calculated crystal stability fields. 



from the pressure isochores reported in Ref. [29(, and 
plot the location in the p-T plane in Fig. QTJ We see that 
it falls within the coexistence region of two crystals of 
significantly different p, S and LDT. This is similar to 
the case of, e.g., the TIP4P2005 model of water [H HI 
and is consistent with the idea that L-L phase separation 
is possible when there is a strong coupling between energy 
and density 50]. 



IV. DISCUSSION 

We calculate the coexistence temperature (along an 
isobar) or pressure (along an isotherm) of two phases by 
determining the point at which the chemical potential 
of those two phases cross, and estimate the uncertainty 
by accounting for the numerical error, typically arising 
from an integration, in evaluating the various terms in, 
e.g., Eq . P21 IT21 and [DH The errors mostly result in a con- 
stant shift in the curves that, given the small difference 
in slopes of p(P) or p(T) between the liquid and crys- 
tal phases, can lead to a large uncertainty in the cross- 
ing. As a check, after calculating the coexistence curve 
through Gibbs-Duhem integration, for the L-S case for 
example, we determine two additional chemical poten- 
tial crossings along different thermodynamic paths and 
the results show good consistency with the Gibbs-Duhem 
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curve. Another indicator of the quality of the results are 
the degree to which coexistence lines cross at the L-S- 
HDT and L-S-LDT triple points. 

Having said this, shifts in the /x(P) or fi(T) curves do 
not affect the slopes, which show in general the first-order 
character of the L-S or L-LDT transitions. For a given 
phase, we determine /i to the point where it is simple 
to determine the equilibrium properties of that phase, 
i.e., to the point where spontaneous transformation does 
not readily occur on the timescale of simulation. We note 
that for the liquid to crystal transitions, the chemical po- 
tential difference between the liquid and crystal at which 
metastability is no longer easily attainable is rather small 
in comparison to other studies |5l[ . Perhaps this is a fea- 
ture of two-dimensional systems, but nonetheless implies 
a very small surface tension if the classical description of 
nucleation is valid. 

The L-LDT and L-S crystallization lines in Ref. [23], as 
noted earlier, were dynamically determined as maximal 
extents of the liquid's ability to exist, and we show here 
that they indeed occur in the metastable liquid. The loss 
of liquid metastability prevents observing any low and 
high density liquids that would exist below the proposed 
L-L critical point because these limit lines radiate from 
the critical point towards higher T. We would like to 
explore the process of crystallization in this vicinity. If 
indeed the L-L critical point proposed for this system 
is obscured by nucleation induced through critical fluc- 
tuations, studying nucleation in the present model may 
help better understanding what may be occurring in wa- 
ter |52l |. 

Notably, for the model at higher P, we provide evi- 
dence for inverse melting, arising from a maximum in P 
in the L-S coexistence line. This phenomenon is rare, and 
seeing evidence for it in such a simple system will allow 
for deeper exploration into the basic physics surrounding 
it. 

The freezing of the liquid to the close-packed solid, i.e., 
the L-HDT transition, appears to be first-order for all T 
that we have explored. For our system size, the free en- 
ergy barrier between the L and HDT basins with p as 
the order parameter at fc^T/e = 5.0 and P = 50.0e/a 2 
is just above lfegT. In the high T limit when the sys- 
tem should behave as hard disks, Mak [53[ has provided 
evidence that the transition should be also first-order. 
Lowering T, the barrier grows and reaches a value of 



~ 2AksT at a simulation conducted on our coexistence 
line at k B T/e = 0.50 and P<r 2 /e = 9.5649 with N = 986, 
thus becoming more strongly first-order. An investiga- 
tion into the region of HDT close to melting would be 
warranted, as Mak has shown that for hard-disks, freez- 
ing occurs to a crystal in which an orientational order 
parameters scales with system size in a way consistent 
with the hexatic phase. 

V. CONCLUSIONS 

We compute a phase diagram using various free en- 
ergy techniques of a two-dimensional SSSW model that 
has been previously shown to exhibit liquid-state anoma- 
lies often associated with the presence of a metastable 
L-L critical point [l8( . We find two low-T crystal phases 
not previously reported. All transitions, including melt- 
ing lines, appear to be first-order for our system size 
of ~ 1000 particles. Thus, it appears that the liquid 
anomalies present in the system do not arise as a re- 
sult of quasi-continuous freezing, as has been previously 
suggested [54|. Previously reported crystallization lines 
fall within respective phase stability regions reported 
here. Interestingly, the difference in chemical potential 
between liquid and crystal phases at the limit where the 
metastable phase can be readily observed is rather small, 
/3A^ - 0.01. 

The L-S coexistence curve exhibits both a maximum 
temperature, indicating that at higher pressure the crys- 
tal is less dense than the melt, and a pressure maximum, 
which means that inverse melting should occur in a spe- 
cific pressure range. Given the scarcity of systems ex- 
hibiting inverse melting, the present model presents the 
opportunity to study this rare phenomenon in more de- 
tail. 
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